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ABSTRACT 

The  following  model  problem  examines  the  thcrmocapillary  feedback  mechanism  im¬ 
portant  at  the  edge  of  weld  pools  and  other  materi-  processes.  A  pool  of  liquid  with  a  flat 
horizontal  free  surface  is  bounded  on  one  side  In  a  vertical  solid  wall,  which  is  maintained 
at  a  cold  temperature  to  unit  depth,  and  at  a  warmer  temperature  below;  far  away  the 
fluid  is  at  the  warmer  temperature.  Surface  tension  is  a  decreasing  function  of  tempera¬ 
ture,  so  that  the  surface  thermal  gradient  drives  flow  toward  the  corner.  W'hen  convection 
is  vigoi'ous,  the  flow  compresses  the  thermal  gradient  which  is  driving  the  flow:  this  posi¬ 
tive  feedback  results  in  small  local  length  scales  and  high  velocities  near  the  corn<^r.  This 
problem  is  examined  through  a  detailed  scaling  analysis  and  through  numerical  s.mulation 
for  a  range  of  parameters. 
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THERMOCAPILLARY  FLOW  NEAR  A  COLD  WALL 


1.  INTRODUCTION 

In  the  processing  of  materials,  often  material  is  melted  and  resolidified.  Several  practi¬ 
cal  processes,  e.g..  welding,  float-zone  purification,  and  Czocliralski  crystal  growth,  involve 
a  pool  of  molten  metal  with  a  free  surface,  with  strong  temperature  gradients  along  the 
surface.  Convection  in  the  molten  metal  is  typically  vigorous  and  significant  to  the  results 
of  the  process,  in  that  it  affects  the  size  and  shape  of  the  pool,  the  heat  transfer,  the 
mixing  of  solutes,  and  ultimately  the  microstructure  of  the  finished  product.  The  forces 
driving  the  convection  include  the  variation  of  surface  tension  with  temperature  along 
the  surface  {thermocapillary  forces),  buoyancy  forces  due  to  thermal  (and/or  solutal)  ex¬ 
pansion,  and  electromagnetic  forces  in  the  case  of  arc  welding  or  electron  beam  welding. 
However,  in  many  cases  (e.g.  laser  welding)  thermocapillary  forces  ])rcdomiuatc.  and  even 
in  cases  where  other  forces  are  stronger  overall,  there  arc  still  important  regions  where  the 
thermocapillary  forces  are  dominant  (i.{'..  cold  corner  regions:  see  Chen.  1987). 

Consequently,  there  have  be<‘n  many  theoretical  studies  of  thermocapillary  flows.  j>ii- 
marily  numerical,  and  a  few  analytical  (review's  arc  given  by  Ostrach.  1982.  and  Davis. 
1987).  Cowley  and  Davis  (1983)  analyzed  the  (tw’o-dimcnsional)  thennocapillaiy  flow'  near 
a  hot  wall  for  vigorous  flow'  (large  Marangoni  number):  here  the  fluid  flows  iq)  the  wall 
then  turns  and  flows  aw'ay  along  the  free  surface,  so  this  w'o\dd  be  called  the  hot  corner 
problem.  The  numerical  studies  of  Zebili  ct.  al.  (1985)  of  flow’  in  a  rectangular  pool  (2-D) 
w’ith  one  hot  and  one  cold  wall,  how'cver.  show'  that  for  moderate  to  small  Prandtl  numbers 
(e.g.,  metals)  the  cold  corner  region  has  by  far  the  strongest  effect  on  both  the  flow'  and  the 
heat  transfer.  This  result  gives  a  different  overall  scaling  than  that  of  Cow’ley  and  Davis, 
although  their  local  hot-corner  scaling  was  validated.  Other  numerical  studies  (e.g..  Zchr 
et.  al..  1987).  when  a  sufficiently  fine  mesh  is  used,  show'  similar  strong  flow  at  the  cold 
corner,  where  the  flow  along  the  free  surface  toward  the  cold  w'all  compresses  the  thermal 
gradient,  thereby  enhancing  the  flow  and  the  heat  transfer.  Great  care  is  nc<cssary  to 
insure  that  the  small  length  scales  of  this  corner  region  are  resolved  numerically:  this  is 
not  always  the  case  (as  noted  by  Chen.  1987). 

Therefore,  it  is  imperative  to  develop  a  theoretical  understanding  of  the  dynamics 
of  the  cold  corner  region.  Being  a  region  of  intense  heat  transfer,  the  details  of  the  flow’ 
can  affect  the  shape  of  the  melt  pool  and  the  cooling  rate,  thus  the  microstructure,  of 
the  material.  At  the  least,  the  dependence  of  the  length,  velocity  and  thermal  scales 
on  the  parameters  (Marangoni  number.  Prandtl  number,  Capillary  number)  needs  to  be 
understood  in  order  for  realistic  numerical  models  to  be  designed  in  a  way  to  resolve  the 
details  in  this  important  region.  But  as  yet.  such  understanding  is  lacking.  In  fact,  in  a 
recent  review,  M.  M.  Chen  (1987.  p.552)  states,  “It  would  seem  then  that  the  stiucture  of 
the  cold  corner  flow  is  one  of  the  most  critical  issues  to  be  studied  in  the  future." 

To  analyze  the  behavior  of  the  cold  corner  region  without  all  the  complications  of  the 
complex  geometry,  phase  change,  and  time  dependence  inherent  in  real  materials  process¬ 
ing  applications,  a  simplified  model  prol)lem  w'ill  be  considered,  much  like  that  of  Cow'ley 
and  Davis  (1983),  as  follows.  A  pool  of  liquid  has  a  horizontal  free  surface  ending  at  a 
vertical  wall,  and  the  upper  .section  of  the  wall  is  cooled:  the  resulting  thermal  gradient 
drives  thermocapillary  flow  towards  the  cold  corner.  The  depth  and  width  of  the  pool  are 
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assuuied  large  compared  to  all  local  length  scales  (which  is  reasonable  for  practical  situa¬ 
tions  with  high  Marangoni  numbers),  so  the  pool  appears  semi-infinite  both  horizontally 
and  vertically. 

This  simplified  problem  is  still  complicated,  and  contains  most  of  the  filatures  of  the 
cold  corner  regions  in  practical  processes,  e.g..  welds.  The  mi.‘;.sing  features  art-  phase 
change  and  surface  deflection,  both  of  which  could  modify  the  geometry  locally  (curved 
wall  and  surface;),  but  are  unlikely  to  change  the  basic  structure  and  scaling. 

The  structure  of  the  thermal  and  flow  fields  will  be  examined  both  through  a  detailed 
scaling  analysis  to  determine  the  dependence  on  the  ])arameters.  and  through  numerical 
simulation  using  two  different  methods. 

2.  PROBLEM  STATMENT 

A  pool  of  incompressible  Newtonian  fluid  is  bounded  on  the  left  by  a  vertical  solid 
w'all.  with  the  upper  portion  of  the  wall  (to  depth  d)  maintained  at  a  cold  temi)erature  Tc- 
while  the  rest  of  the  wall  is  at  the  hot  ambient  temperature  T/,  of  undisturbed  fluid  far  from 
the  corner.  (Sec  figure  1.)  Above  the  horizontal  free  surface  of  the  liquid  is  an  inviscid. 
nonconducting  gas.  Surface  tension  is  assumed  strong  enough  to  keep  the  free  surface  flat 
(small  Capillary  number),  but  with  surfa<-e  tension  variations  due  to  a  linear  d('i)endonce 
on  temperature.  The  resulting  flow  is  assumed  to  be  two-dimensional  and  steady. 

Then  the  equations  governing  the  thcrmocapillary  convection  in  the  cohl  corner  are 
conservation  of  miuss.  momentum,  and  energy: 

V-u  =  0  (2.1) 

pu  ■  V  u  = —V p  + /I  V‘ u  (2.2) 


pCpXX 

■  VT 

=  kV^T 

(2.3) 

with  the  boundary  conditions: 

at  (/  =  0  : 

T  = 

0.  V  =  0.  pity  =  qTx 

(2.4a.b.c) 

at  X  =  0  : 

T  = 

V  A 
a.  a. 

11 

II 

(2.5a.b.c) 

as  x.  y  —*  oc  : 

T  ^ 

T/,.  u.  V  — ►  0 

(2.Ca.b.c) 

Here  u  is  the  velocity  vector  with  components  u  and  v  in  the  x  (horizontally  rightward) 
and  y  (vertically  downward)  directions,  p  is  prc.s.surc.  T  is  tcmi)crature.  p  is  density,  p  is 
viscosity.  Cp  is  specific  heat,  k  is  thermal  conductivity,  and  7  (assumed  constant  and  neg¬ 
ative)  is  the  derivative  of  the  surface  tension  with  resi)ect  to  temi)erature.  The  boundary 
conditions  specify  that  the  wall  is  piecewise  isothermal  with  no  fluid  slip,  and  the  flat  free 
surface  is  thermally  insulated,  with  thcrmocapillary  forcing. 

The  equations  can  be  nondimensionalized  by  scaling  lengths  by  d.  temperature  dif¬ 
ferences  by  AT  =  Th  —  Tc.  and  velocities  by  m,,  =  7  AT//c.  The  resulting  dimensionh'ss 
equations  are: 

V  •  u  =  0 

i?,u  •  V  u  =  —V p  -f  u 
Mu- VT  =  V'T 
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(2.7) 

(2.8) 
(2.9) 


with  the  boundary  conditions: 


at  y  =  0  : 

Ty  =  0.  V  =  0.  Uy  =  Tj 

(2.1()a.b.c) 

(  -1-  V  <  1 

at  X  =  0  : 

T= \  •  .  v= V = 0 

(2.11a.b.c) 

1  0.  y>l 

X.  y  —>  cc  : 

T  — +  0.  ti.  t'  — »  0 

(2.12a.b.c) 

where  u,  etc.  from  here  on  denote  the  dimensionless  quantities.  The  two  dimensionless 
parameters  are  the  Marangoni  number  M  =  Uad/K  and  the  Reynolds  number  R  =  Uadfy. 
where  k  =  kj pcp  is  the  thermal  diffusivity  and  u  =  pfp  is  the  kinematic  viscosity.  Their 
ratio  gives  the  Prandtl  numix'r;  P  =  v [n  =  M/R. 

For  the  numerical  solutions  below,  it  is  convenient  to  eliminate  the  pressure  by  ado])t¬ 


ing  a  stream-function/vorticity  formulation  for  the  flow: 

/?  u  •  V  =  V*  u.’  (2.13) 

(2.14) 

V  =  '^y  .  V  =  — (2.1oa.b) 

where  is  the  stream  function  and  cj  is  the  vorticity.  with  the  flow  boundary  conditions 

at  y  —  (]  :  <!»  =  =  0.  a>  =  — Tx  (2.1Ca.b.c) 

at  X  =  0  :  = 'J'y  =  0  (2.17a.b.c) 

as  x,y —*  oc  :  (2.18a.b) 


3.  SCALING  ANALYSIS  AND  REGIMES  OF  BEHAVIOR 

The  structure  of  the  thermal  and  flow  fields  can  take  on  different  forms.  dei>ending  on 
the  values  of  the  two  governing  j)aramctcrs.  the  Marangoni  number  M.  wdiich  measures  the 
importance  of  thermal  convection  relative  to  thermal  diffusion,  and  the  Prandtl  number 
P.  which  is  the  ratio  of  viscous  to  thermal  diffusion,  a  material  pro]>erty.  (Equivalently, 
one  could  use  Reynolds  number  R  =  M/P  as  the  second  jjarameter.  the  ratio  of  iiun  tial  to 
viscous  forces.)  Here  we  derive  the  appropriate  dimensionless  scales  for  the  four  different 
asymptotic  regimes  of  behavior. 

To  examine  the  dominant  balances  in  the  cohl  corner,  three  scales  suffice:  the  hori¬ 
zontal  length  scale  I  for  the  thermal  gradient  along  the  surface,  the  vertical  viscous  length 
scale  6  for  the  velocity  shear  at  the  surface  (which  turns  out  to  be  the  same  scale  as  for 
the  velocity  shear  on  the  wall),  and  the  velocity  scale  U  for  flow  along  the  surface.  The 
vertical  thermal  length  scale  is  determined  by  the  boundary  condition  on  the  wall,  and  so 
is  0(1).  The  thermocapilh'iry  stress  condition  (2.10c)  scales  as: 


so  that  U  ~  6/1.  In  the  energy  and  vorticity  equations  (2.9  and  2.13).  the  terms  for 
convection  in  each  direction  scale  the  same,  but  not  so  for  diffusion,  and  by  (2.1G).  the 
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surface  vorticity  scale  is  1//.  so  (2.9)  and  (2.13)  scale  as: 


U 


P 

1 


+  1 


—  /«>w  —  -4-  — ' 

p  P 


(3.2) 

(3.3) 


For  small  enough  M,  thermal  convection  is  negligible,  implying  /  ~  1.  and  the  thermal 
field  is  essentially  conductive,  decoupled  from  the  flow.  But  for  large  enough  M  convection 
becomes  important,  and  the  strong  surface  flow  toward  the  wall  compresses  the  thermal 
gradient  along  the  surface,  which  in  turn  strengthens  the  local  driving  force  for  the  flow. 
This  reduces  the  horizontal  thermal  length  scale  I  to  the  point  that  thermal  diflusion  away 
from  the  wall  balances  conv('ction  toward  the  wall,  so  the  local  cffrrUvr  Maraugf)ni  numlx'i 
is  order  unity:  Me//  =  MUl  ~  1.  Then  the  externally  imposed  length  scale  (dimensional 
d  above)  is  no  longer  directly  relevant  to  the  com])ressed  cold  corner  region.  (In  this  case, 
the  local  importance  of  inertia  is  better  indicated  by  whether  viscous  or  tln'iinal  diflusion 
is  more  efficient,  i.e..  by  P  rather  than  R.) 

Similarly,  for  small  enough  R.  inertia  is  negligible  everywhere,  implying  (S  ~  /.  and  the 
flow  is  dominated  bj'^  viscous  forces.  For  large  enough  R  inertial  forces  beconu'  dominant 
and  viscous  effects  are  confined  to  boundary  layers  of  thickness  6  I  along  the  surface  and 
the  wall,  where  the  local  effective  Reynolds  number  Rcff  =  RUb~  jl  ~  1.  (Both  layers  are 
of  comparaltle  thickness  because  the  pressiire  field  outside  the  layers  has  the  same  length 
scale  in  both  directions.)  The  above  gives  the  scaling  for  each  regime. 

When  the  thermal  field  is  conductive  and  the  flow  field  dominat  ed  by  viscous  forces 
[M  -C  1  and  R  1.  or  P  >  ilf).  all  three  scales  are  of  order  unity:  I  ~  1.  ~  1.  U  ~  1. 
Thus  in  this  case  (only),  the  scaling  \i.sed  in  the  nondinn'iisionalization  is  a])])i()priat(' 
everywhere.  Within  this  regime,  the  solution  is  fully  two-dimensional  with  no  fine  structure 
and  is  nearly  independent  of  the  parameters. 

For  the  conductive  case  with  inertial  flow,  the  additional  resistance  of  inertia  reduces 
both  the  velocity  scale  and  the  viscous  length  scale:  C7  ~  6  ~  R~ (while  /  ~  1  still). 
This  reduced  velocity  also  reduces  the  effective  Marangoni  number,  such  that  this  regime 
applies  when  M  RK  or  M  •C  P~K  with  i?  ^  1,  or  P  <C  M.  (Note  that  this  gives 
the  same  boundary  layer  scaling  as  Zebib  et  al..  1985,  cxcci)t  they  made  an  error  on  the 
passive  wall  layers,  as  pointed  out  by  Chen.  1987.)  Here  the  vorticity  generated  by  the 
shear  stresses  on  the  surface  and  the  wall  are  confined  to  the  thin  boundary  layers. 

When  thermal  convection  is  imi)ortant  but  inertia  is  not  (M  >>  1  and  P  1).  surface 
thermal  varitions  are  compressed  to  a  narrow  region,  beyond  which  the  thcrmocainllary 
forcing  is  small,  so  ^  ~  ~  M~^  and  ~  1.  However,  the  strong  inward  flow  along 

the  surface  turns  downward  and  away  from  the  no-slip  wall  (and  weakens  rapidly  with 
distance),  such  that  no  thermal  boundary  layer  is  formed  on  the  wall:  rather,  vertical  and 
horizontal  variations  arc  comparable. 

The  most  important  regime  for  materials  processing  is  where  thermal  convection  is 
important  and  P  1.  the  latter  being  generally  true  for  metals.  In  this  case,  within  the 
compressed  thermal  region  there  are  thin  viscous  boundary  layers  along  the  surface  and 
wall.  Then  I  ~  M~^ P~K  d  ~  M“^.  and  U  ~  P*.  i.e.,  the  additional  resistance  of  inertia 
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decrejises  the  velocity  scale  and  thus  iiicieiises  the  thermal  scale  by  a  factor  of  \/  P  relative 
to  the  purely  viscous  case.  Again  the  reduced  velocity  changes  the  thermal  convection 
scaling,  and  large  Marangoni  number  here  means  Af  ^  P~^ . 

The  approximate  divisions  between  the  four  asymptotic  regimes  are  shown  in  figure 
2.  For  a  material  of  small  P.  as  M  is  increased  from  zero,  at  first  the  temperature  field  is 
conductive  and  the  flow  dominated  by  viscous  forces,  then  the  flow  becomes  primarily  iner¬ 
tial  and  viscous  boundary  layers  form,  and  finallj'  thermal  convection  becomes  dominant, 
shrinking  all  local  length  scales  in  the  corner. 

4.  VISCOUS  CORNER  REGION 

There  is  a  region  in  the  corner,  for  any  A/  and  P.  where  viscous  stresses  from  the  wall 
limit  the  flow  and  both  inertia  and  thermal  convection  are  negligible,  so  the  temperature 
is  a  linear  function  along  the  surface.  Locally  the  thermocapillary  stress  is  constant,  and 
the  flow  is  given  by  a  similarity  solution  (Moffatt.  19C4.  although  the  published  version  is 
incorrect). 

If  the  flat  free  surface  makes  an  angle  o  with  the  solid  wall,  then  a  constant  unit 
surface  shear  stress  toward  tin?  corner  giv(‘s 

T,  [{sin  2o  —  2o  )(cos2^^  —  1 )  —  (cos 2o  —  1  )(sin  2^  —  2^)1 

=  \ ^ ^  (4.1) 

4  sill  2a  —  2a  cos  2a 

where  r  and  6  are  polar  coordinates,  with  0  increasing  from  the  wall  to  the  free  surface 
(see  figure  1).  In  the  special  case  here,  a  =  and 


^(7'.  0)  =  r* 

u  =  r  I 


-(1  -cos2^)-l  — (sin2^-20) 
4  27r 


-  sin  20  -I — (cos 26*  —  1) 

2  TT 


^(1  -  cos2t?)  -f  -(sin  2^  -  20) 

2  TT 


a;  =  1 - 0 

IT 

4  , 

p-po  =  -hir 
TT 


(4.2a) 


(4.2b) 


(4.2c) 


(4.2d) 


where  r  and  0  are  unit  vectors  in  the  coordinate  directions,  and  Pq  is  some  reference 
pressure.  Figure  3  shows  the  streamlines  for  this  similarity  flow.  The  zcro-vorticity  contour 
extends  from  the  corner  at  a  45°  angle,  dividing  the  negative  vorticity  on  the  surface  from 
the  positive  vorticity  on  the  wall. 

This  is  the  form  of  the  flow  in  the  cold  corner  on  the  smallest  length  scale,  where 
all  the  above  flow  quantities  would  be  multiplied  by  ri(O.O)  (which  scales  as  1//).  The 
velocity  grows  linearly  with  distance  r  from  the  corner  and  the  local  length  scale  is  r.  which 
can  be  used  to  estimate  the  range  over  which  the  similarity  solution  is  applicable.  For  the 
two  viscous-flow  regimes  mentioned  above,  the  linear-temperature  approximation  requires 
r  •C  /,  so  for  the  conductive  regime  r  -C  1  while  for  the  convective  regime  r  <C  AI~^. 
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For  the  two  inertial-flow  regimes,  the  more  restrictive  condition  is  that  locally  inertia  is 
negligible,  and  the  local  velocity  scale  is  r/l.  so  we  must  have  ///?:  in  the  conductive 

regime  this  implies  r  <C  .  while  in  the  convective  regime  tlu-  range  is  r  <  or 

r//  < 

This  gives  an  estimate  of  the  resolution  retpiired  for  a  numerical  model  to  resolve  all 
the  details  of  the  cold  corner  flow.  Because  of  the  corner  singularity  in  u  and  p.  sirectral 
methods  would  be  inappropriate:  instead,  finite-difference  or  finite-element  methods  coukl 
be  used.  Then  if  the  first  grid  point  is  in  the  similarity  range,  no  details  of  the  velocity  field 
in  the  corner  will  be  lost.  In  addition,  the  similarity  form  may  be  useful  as  a  ■■matching " 
type  boundary  condition  for  the  singularity  at  the  origin. 

5.  GREEN'S  FUNCTION  METHOD  FOR  VISCOUS  CASE:  /?  --  0 

When  the  Reynolds  number  is  sufficiently  small,  inertia  is  negligible  throughout  the 
flow  field,  and  to  lowest  approximation  the  flow  is  governed  by  the  biharnionic  ecpiation: 

V‘(V2^)  =  0  (5.11 


Then  the  flow  everywhere  depends  only  on  the  instantaneous  thermal  gradient  along  tlu' 
surface  (even  if  the  flow  is  unsteady).  This  allows  the  flow  field  to  be  rei)resented  using 
the  Green's  function  for  a  point  force  near  a  rigid  wall,  directed  toward  the  wall  (see.  e.g.. 
Heusimoto  (Ic  Sano.  1980.  or  Blake.  1971). 

The  Green's  function  for  the  stream  function  due  to  a  unit  (dimensionh’ss)  force  in 
the  negative  x  direction  applied  on  the  surface  at  the  point  (Cl^) 


G{x.y.i) 


(5.2a) 


where  _ 

=  y/{x  +  iy-+y\  r_  =  v'(.r-02  +  r  (5.2b) 

The  flow  field  due  to  such  a  point  force  at  ^  =  1  is  shown  in  figure  4.  (Note  that 
all  distances  scale  with  (,.)  The  far  field  {7-  ^  can  be  used  as  an  approximation  to  the 
viscous  thermocapillary  flow  at  large  enough  distances  that  the  distributed  forcing  of  the 
thermal  gradient  can  be  replaced  by  a  j>oint  force  near  the  wall.  It  can  be  shown  that  in 
the  far-field  the  zero  vorticity  contour  approaches  an  angle  of  ^  from  the  surfact'.  As  r 
increases,  the  velocity  decays  rapidly  as  u  ~ 

This  rapid  decay  of  velocity  has  the  consequence  that  no  thermal  boundary  layer 
forms  on  the  wall.  While  this  result  only  applies  directly  when  inertia  is  negligible,  one 
would  expect  that  the  same  order  of  decay  applies  to  the  inertial  case  in  the  far  field,  by 
analogy  with  the  Squire-Landau  jet  (Landau,  1944.  and  Squire,  1951)  where  the  velocity 
due  to  an  isolated  point  force  decays  like  u  ~  0(r~^)  for  all  Reynolds  numbers.  Note  also 
that  near  the  wall  (except  near  the  surface),  the  flow  is  away  from  the  wall.  Hence,  though 
the  surface  temperature  gradient  may  be  highly  compressed  from  vigorous  convection, 
there  is  no  thermal  boundary  layer  along  the  wall,  i.e.,  vertical  and  horizontal  temijerature 
variations  are  expected  to  be  of  the  same  order.  This  is  in  marked  contrast  with  the 
seemingly  comparable  case  of  a  hot  corner.  There  the  flow  tends  to  exi)and  the  surface 
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tlu'i'iual  gradient  over  long  distaiu'c.  so  the  forcing  is  similar  to  a  distrihutial  liiu'  force, 
and  thermal  boundary  layers  foru)  both  on  the  wall  and  the  surface  (Cowley  and  Davis. 
1983).  This  difference  between  concentrated  forcing  and  distribut«!d  forcing  changes  tin 
fundamental  character  of  tlu'  flow. 

Using  the  Green's  function,  tlnui  the  stream  function,  velocity  components,  and  voi- 
ticity  in  the  viscous  thcrmoc.apillary  flow  are  given  by: 


'J'(x.y) 
uix.y) 
vix.  y) 
uj{  x.  y  ) 


i  /  T'iO 

/» 

1 


2.r<l 

l-'J 

Jo 

—  f  T' 

2" 


iu(  j;±)+  r 


,.2 


1 


2./-^ 


1  - 


■T  +  ^  X-(, 


(O  2// 


2y^ 

o 

m-r  +  Oy 

■'i 


1  - 


2.T(.r  +  0 
7-i 


( 5. 3a ) 
(l^  (o.3b) 
(o.3e) 
(o.3d ) 


using  the  shortened  notation  T'(.r)  =  T,  (.7-.9). 

To  calcidate  the  steady  flow  for  various  values  of  Marangoni  number  (with  7?  =  0).  the 
time-dependent  equation  for  the  tcmi)eraturc  field  was  exidicitly  intc'grated  in  time  until 
the  steady  state  was  reached.  A  finite  sejuare  domain  on  a  uniform  grid  was  used,  where  on 
the  artificial  boundaries,  diffusion  across  the  boimdaries  is  neglected  and  where  convection 
is  inward  the  fluid  outside  is  assumed  to  Ix'  at  zero  temperature.  The  effects  of  these 
artificial  boundaries  on  the  cold  corner  ari'  presumably  small  if  the  l)oundari(\s  ;u'e  sc'veral 
units  away,  since  the  thermal  field  decays  quickly  with  distance.  Central  finite'  difie're-ncc's 
were  used,  with  upwind  differencing  for  the  convective  terms.  The  velocity  at  each  j)oint 
was  evaluated  from  the  above  inte'grals  {5.3b.c).  using  an  analytic  aj)j)roximation  arceiind 
the  singularity  at  (.r.  0).  and  elsewhere  using  the  trapezoidal  rule  with  first-order  differi'iict's 
for  the  thermal  gradient,  as  central  differences  here  resulted  in  a  numerical  instability. 

The  residts  for  a  range  of  M  are  shown  in  figure  5.  The  transition  from  the  i)rimarily 
conductive  regime  to  the  convective  regime  as  M  increases  is  a])part'nt.  As  the  surface 
gradient  becomes  more  compact,  the  strong  downward  convection  away  from  the  corner 
extends  the  thermal  region  into  a  sort  of  broad  diagonal  plume;  no  thermal  boundary  layer 
forms  even  for  large  M.  (Where  the  isotherms  intersect  the  artificial  boundaries,  slight 
local  distortion  due  to  the  artificial  boundary  conditions  is  ai)parent.)  For  large  M.  the 
flow  outside  the  cold  corner  is  qualitatively  similar  to  the  flow  from  a  jroint  source  (sc(' 
figure  4). 


G.  NUMERICAL  SOLUTIONS  FOR  FULL  SYSTEM 

When  inertia  is  not  negligilde.  the  full  coupled  nonlinear  equations  must  be  .solved. 
Unlike  the  Green  s  function  approach,  this  requires  flow  boundary  conditions  at  the  arti¬ 
ficial  boundaries  of  the  computational  domain:  here  these  boundaries  are  assumed  to  be 
impermeable  and  shear-free,  as  well  as  isothermal.  These  conditions  constrain  both  the 
thermal  and  flow  fields  (compared  to  the  Green's  function  method),  enhancing  recircula¬ 
tion  and  preventing  long  thermal  plumes.  Still,  witli  the  artificial  bomidarics  several  units 
away,  their  effects  on  the  cold  corner  are  expected  to  be  small. 


7 


AgaiiK  the  numerical  method  involves  stepping  the  unsteady  ecjuations  forward  in  time 
until  steady  state  is  reached.  At  each  time  stej).  the  convection-diffusion  etpiations  for  tem¬ 
perature  (2.9)  and  vorticity  (2.13)  are  .solved  by  the  Alternating-Direction  Imjdieit  (ADI) 
scheme,  where  the  convective  terms  are  evaluated  by  the  Eulerian-Lagrangiaii  Method 
(ELM.  see  Cheng  et  al..  1984)  using  the  velocity  field  from  the  old  stream  function  ainl 
uiwind  bilinear  interpolation.  (The  ADI  method  avoids  diffusive  numerical  instability,  so 
the  time  step  is  only  limited  by  convective  stability.)  After  several  tinu'  st(*ps.  th*'  Poisson 
equation  for  the  stream  function  (2.14)  is  solved  by  Gauss-Seidel  iteration  with  Succes¬ 
sive  Over- Relaxation  (SOR).  Steady  state  is  assumed  when  the  pointwise  RMS  change  in 
stream  function  is  below  a  certain  tolerance. 

A  non-uniform  Cartesian  grid  wjus  employed.  (The  program  allows  arl)itrary  s])a(  in'.; 
of  points  in  each  direction.)  Using  a  three-point  difference  scheme,  only  first-order  a«  (  urac\ 
is  possible  for  the  second  derivatives:  the  differencing  em])loycd  becomes  second-onh'r  in 
the  limit  of  uniform  spacing.  Specifically,  the  following  one-dimensional  difference  formulae 
wer(‘  us('d  (derivc'd  from  Taylor  series): 


4-  d:r_)  ^  dx^dx^  ^  dx^{dx^  4-  r/.v_) 
-f  0{f'"dx+dx.) 

2-2  2 
dx^{dx^  +  dx-)  dx+dx-.  ^  dx+{dx+  4-  d.7_) 

4-0(r'(r/:r+  -d:r_)) 


(C.la) 


(r-.ii'i 


where  r/.T_  and  dx+  are  the  distanci's  to  the  grid  points  below  and  ahovc'  the  cuna'iit  point, 
with  /_  and  the  corresponding  values  of  the  function. 

The  grid  spacing  in  each  direction  was  chosen  to  have  a  clos(>ly  i)ack('d  region  of 
uniform  spacing  by  the  surface  or  wall,  a  widely  s})aced  region  of  uniform  s])aring  out  near 
the  artificial  boundaries,  and  in  betw'ccn  a  region  of  smoothly  (exi)onentially)  changing 
spacing.  This  was  generated  by  apjdying  the  following  function  to  a  uniform  grid  (in  C 
say): 


X 


=  f(0  = 


In  v^jD.  <  s  <  si 

1  -  r  Inr  (1  -  ^)/D.  ^2  <  ^  <  1 


(G.2) 


where  r  is  a  parameter  giving  the  ratio  of  outer  to  inner  sjiacing.  and  ^2  jiarametc'is 
delineating  the  three  regions.  D  =  {r  —  1)(^2  —  )  +  lnr(7-  4-  ^vxu  1^^  fl>(‘ 

position  of  the  artificial  boundary. 

It  was  found  that,  even  when  the  time  stej)  easily  satisfied  convective  stability  require¬ 
ments.  nonetheless  instabilities  developed  in  the  vorticity  near  the  wall.  Several  diffcrcMil 
formulations  for  the  wall  boundary  condition  on  vorticity  were  tried,  to  no  avail.  However, 
by  under-relaxing  the  changes  in  wall  vorticity  only  during  the  initial  adjustme  nt  jieriod. 
the  instability  was  eliminated.  The  formula  used  to  calculate  the  vorticity  at  the'  wall 
(without  under-relaxation)  from  the  stream  function  is 


UJf) 


-  ^2) 

1-2  _  1-2 
•*'2  •'1 


+  0('I'"'a-2) 


(G.3) 


8 


where  subscripts  0,  1,  and  2  refer  to  the  wall  and  the  first  two  grid  i)oiiits.  While  only 
first-order,  this  forniulation  is  independent  of  avoids  re-using  the  boundary 

condition  for 

'  Results  are  shown  in  figures  6  and  7  for  P  =  1  and  P  =  0.01.  for  a  variety  of  A/.  The 

conductive,  viscous  regiinc  is  represented  by  the  case  A/  =  0.01.  P  =  \  (figure  Ca).  For 
P  =  1.  with  increasing  M.  the  surface  thermal  gradient  becomes  compressed  (similar  to 
the  P  =  DC  case  computed  by  Green's  functions)  along  with  the  vorticity  on  the  surfac  e. 
However,  inertia  is  no  longer  negligible,  and  so  the  flow  down  the  wall  has  less  of  an 
outward  component.  Also,  the  artificial  boundaries  modify  the  ‘‘plume."  keeping  it  from 
the  boundaries  and  turning  it  ujnvard  due  to  the  recirculation.  For  P  =  0.01  (as  typical  in 
liquid  metals),  with  increasing  M  inertial  effects  confine  the  vorticity.  forming  clear  viscous 
boundary  layers,  before  thermal  convection  becomes  strong.  Also,  a  counter-rotating  cell 
forms  in  the  lower  part  of  the  domain.  As  thermal  convection  compresses  the  surface 
gradient,  tlie  surface  viscous  boundary  layer  remains  limited  to  the  cold  corner.  Note  that 
for  large  M  and  small  P,  the  thermal  length  scale  is  small  compared  to  the  overall  domain, 
and  the'  viscous  length  .scale'  is  small  com[)ared  to  the  thermal  scale.  im])osing  severe  (local) 
resolution  requirements  on  any  numerical  model. 

The  numerical  results  for  P  =  1  and  P  =  0.01  are  compared  with  the  scaling  analysis 
in  figure  8.  The  thermal  length  scale  (or  rather  1//)  is  estimated  from  the  thermal  gradient 
Tr  at  the  wall  (based  on  the  first  grid  point  along  the  surface).  The  velocity  scale-  {U) 
is  estimated  from  the  maximum  velocity  Umar  at  a  grid  point  on  the  surface  (though  the 
actual  maximum  might  be  expected  to  fall  between  grid  points).  The  viscous  length  scale 
(A)  is  approximated  by  the  position  Xr,utx  of  the  grid  point  witJi  the  maximum  velocity. 
For  P  =  1,  the  transition  is  apparent  from  the  conductive  viscous  regime  (where  L  U. 
and  S  are  nearly  independent  of  M)  to  the  convective  viscous  regime  (where  /  dc  A/~^  and 
6  DC  A/"\  with  U  roughly  constant).  For  P  =  0.01,  two  transitions  are  seen,  from  the 
conductive  viscous  regime  to  the  conductive  inertial  regime  {I  nearly  constant.  U  a  . 

and  6  oc  to  the  convective  inertial  regime  (/  dc  M~^ .  b  dc  M~^  .  U  nearly  constant ). 

7.  DISCUSSION  AND  CONCLUSIONS 

The  practical  importance  of  thcrmocapillary  convection  in  materials  processing,  along 
witli  the  complications  inherent  in  tyihcal  processes  (e.g.  curved  ijitcrface.s.  pliasc  change, 
etc.),  insure  that  numerical  simulations  will  remain  one  of  the  main  theoretical  tools  for 
understanding  such  systems.  This  work  predicts,  a  priori,  the  resolution  requirements  for 
such  numerical  models  to  accurately  represent  the  high  heat  transfer  and  rapid  velocity 
variations  in  the  cold  corner  region. 

The  structure  of  the  corner  depends  on  two  dimensionless  parameters  indicating  the 
driving  force  for  convection  and  relative  importance  of  viscosity:  the  Marangoni  number 
»  M  (based  on  the  overall  temperature  difference,  overall  thermal  length  scale,  and  material 

properties)  and  the  Prandtl  number  P  (a  material  property).  Hence  there  are  four  asymp¬ 
totic  regimes  (shown  in  figure  2)  depending  on  whether  thermal  convection  and  inertial 
*  forces  are  locally  important.  For  large-P  materials  (e.g.,  organics),  the  corner  flow  is  dom¬ 

inated  by  viscous  forces,  and  for  large  M  the  surface  thermal  variations  are  compressed, 
i.e..  the  local  length  scale  is  decreased  (and  heat  transfer  increased).  For  small-P  materi¬ 
als  (e.g..  metallics),  inertia  becomes  important  before  thermal  convection  with  increasing 


9 


M .  and  thin  viscous  boundary  layers  form  within  the  thermal  region.  When  M  is  large 
(compared  to  three  levels  of  length  scales  must  be  resolved  (overall,  thermal,  and 

viscous),  a  severe  requirement  on  numerical  models. 

The  scaling  was  derived  from  a  sinqile  problem  designed  to  isolate  the  feedback  nu'ch- 
anism  of  the  cold  corner.  Numerical  calculations  by  two  methods  for  closi'ly  related  prol>- 
lems  (where  the  corner  is  necessarily  less  isolated  due  to  the  finite  domain)  illustrate  the 
changing  structure  of  the  cold  corner  and  show  the  details  of  the  transitions  between  the 
asymptotic  regimes.  The  numerical  results  are  consistent  with  the  scaling  analysis. 

One  surprising  result,  in  contrast  with  the  hot  corner  problem  of  Cowley  and  Davi^ 
(1983).  is  that  no  thermal  boundary  layers  form.  This  difference  is  due  to  the  forcing  being 
limited  to  a  relatively  concentrated  region  in  the  cold  corner,  while  for  the  hot  corner  th(> 
forcing  is  distributed  over  an  extended  region  (the  horizontal  thermal  variations  being 
extended  by  convection). 

To  compare  with  the  numerical  results  of  Zebib  et  al.  (1985)  for  P  =  1.  note  tlu' 
different  boundary  conditions:  their  domain  had  a  (flat)  free  stirface.  a  hot  and  a  cold 
no-slip  wall,  and  an  insulated  no-sliit  bottom.  Then  even  when  the  cold  corner  was  highly 
compressed,  there  were  still  significant  thermal  variations  along  most  of  the  surface  due  to 
the  hot  wall,  so  the  overall  scaling  is  similar  to  the  conductive  inertial  case  here.  And  aj)- 
parently  even  their  cold  corner  was  modified  by  the  external  bulk  flow,  for  their  inaxinnun 
vorticity  scaled  as  w  cx  whereas  here  <jJ  oc  M .  (The  resolution  reciuirenu'nts  (huived 

here  can  thus  be  considered  conservative.)  Hence  one  important  question  remaining  is 
under  \vhat  conditions  can  the  cold  corner  be  considered  locally  determined  (a.-  it  is  here). 

In  real  materials  processes,  the  surface  is  free  to  deflect,  and  the  i)o.sition  aiifl  shajie  of 
the  solid-liquid  interface  depends  on  the  thermal  field.  These  effects  greatly  comiflicate  the 
problem  geometrically,  yet  the  dominant  dynamic  balance  should  remain  that  comsidered 
here.  Future  work  will  investigate  these  effects  in  the  cold  corner. 
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FIGURE  CAPTIONS 


Figure  1.  Problem  Formulation:  a  liquid  quarter-space  is  bounded  above  by  a  flat  free 
surface  subject  to  thermocapillary  forcing,  and  is  bounded  on  the  left  by  a  rigid  vei  tical 
wall,  at  temperature  to  depth  d  and  at  the  warmer  temperature  T/,  below,  whieh  i.<  al."*) 
the  ambient  temperature  of  the  undisturbed  fluid  far  away. 

Figure  2.  Asymptotic  Regimes:  the  four  different  asymptotic  behaviors  corresj)oiul  to 
whether  the  thermal  field  is  controlled  primarily  by  convection  or  by  conduction,  and 
whether  the  flow  is  dominated  by  viscous  or  inertial  forces. 

Figure  3.  Viscous  Corner  Flow:  streamlines  (solid,  black)  and  vorticity  contours  (dashed, 
light  gray)  for  the  self-similar  solution  valid  where  the  surface  temperature  is  approximately 
linear  and  inertia  is  negligible. 

Figure  4.  Green's  Function:  streamlines  (solid,  black)  and  vorticity  contours  (dasln'd.  light 
gray)  for  the  non-incrtial  flow  dtie  to  a  point  force  on  the  surface  at  (1.0).  directed  toward 
the  wall.  The  flow  is  recirculating  and  decays  i-apidly  with  distance. 

Figure  5.  Numerical  Results  for  /?  — >  0  (P  — >  oc);  velocity  vectors  and  isotherms  (solid, 
dark  gray)  of  numerical  solutions  irsing  the  Green's  function  fonmdation.  fur  a  range  of 
Marangoni  numbers.  On  the  artificial  boundaries  (bottom  and  right),  normal  thermal 
diffusion  is  neglected  and  incoming  fluid  is  assumed  Lsothermal.  (a)  M  =  10:  (b)  M  =  30: 
(c)  M  =  100;  (d)  M  =  300:  (e)  M  =  1000. 

Figure  C.  Numerical  Results  for  P  =  1:  isotherms  (solid,  dark  gray),  streamlines  (thin, 
black)  and  vorticity  contours  (da.slied.  light  gray),  for  tlie  numerical  solution  of  tin'  full 
coupled  system.  The  artificial  boundaries  (bottom  and  right)  are  assumed  isothermal, 
impenetrable,  and  shear-free,  (a)  M  =  0.01:  (b)  M  =  10:  (c)  M  =  30:  (d)  M  =  100:  (e) 
M  =  300:  (f)  M  =  1000:  (g)  detail  of  M  =  1000. 

Figure  7.  Numerical  Results  for  P  =  0.01:  sec  previous  caption,  (a)  M  =  1:  (b)  M  =  10: 
(c)  M  =  100:  (d)  detail  of  M  =  100:  (o)  M  =  1000:  (f)  detail  of  M  =  1000:  (g)  M  = 
10.000:  (h)  detail  of  M  =  10.000. 

Figure  8.  Summary  of  Numerical  Scales:  wall  temperature  gradient  (diamonds),  maxi¬ 
mum  velocity  Urnax  (triaiiglcs).  and  position  of  maximum  velocity  (stars),  as  functions 
of  Marangoni  number  M.  from  numerical  solutions  of  full  system.  Lines  of  slope  1.  -1.  and 
-1/3  are  shown  for  comparison  with  scaling  analysis,  (a)  P  =  1:  (b)  P  =  0.01. 
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